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We study coherent multiple Andreev reflections in quantum SNS junctions of finite length and 
arbitrary transparency. The presence of superconducting bound states in these junctions gives rise 
to great enhancement of the subgap current. The effect is most pronounced in low-transparency 
junctions, D < 1, and in the interval of applied voltage A/2 < eV < 2A, where the amplitude 
of the current structures is proportional to the first power of the junction transparency D. The 
resonant current structures consist of steps and oscillations of the two-particle current and also of 
multiparticle resonance peaks. The positions of the two-particle current structures have pronounced 
temperature dependence which scales with A(T), while the positions of the multiparticle resonances 
have weak temperature dependence, being mostly determined by the junction geometry. Despite 
the large resonant two-particle current, the excess current at large voltage is small and proportional 
to D 2 . 
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I. INTRODUCTION 



Transport properties of small conducting structures are strongly influenced by size effects. Oscillation of magne- 
toresistance in thin metallic films, and quantization of conductance in narrow wires and point contacts are examples of 
such effects. Size effects in superconducting tunneling have attracted attention since early experiments by Tomaschla. 
In these experiments, oscillations of the tunnel conductance as a function of applied voltage were found for tunneling 
from a superconductor to a thin superconducting film of an NS proximity bilayer. The geometric resonance nature 
of the effect was clearly indicated by the dependence of the period of oscillations on the thickness of the supercon- 
ducting film. Similar conductance oscillations for tunneling into a normal metal film of NS bilayers were reported 
by Rowell and McMillana. Later on an even more pronounced effect - steps on the-current-voltage characteristics 
of SINS junctions at applied subgap voltages, eV < 2A - was observed by Rowellu (for a review see Ref. ||). In 
addition to the dependence on the thickness of the N-film, the period of the current steps also shows temperature 
dependence which scales with the temperature dependence of the superconducting gap A(T). The current steps occur 
at applied subgap voltages, eV < 2A, and they are understood as resonant features due to quasiparticle tunneling 
through superconducting bound states existing in INS wells at energies lying within the superconducting gap, \E\ < A, 
de Gennes — Saint- James levelsO. 

Recently, properties of superconducting bound states have attracted new attention in connection with observations 
of conductance anomalies in mesoscopic NS structures. Resonant-oscillations of the subgap conductance in mesoscopic 
quasi-ballistic NS junctions have been reported by Morpurgo et a!3. JChe zero bias conductance peak, a related resonant 
phenomenon in diffusive NS junctions, was discovered even earlienj'tj and then intensively studied both theoretically 
and experimentally (for a review see Ref. ^|) . 

In these recent studies of mesoscopic junctions, attention was switched from the single-particle current through 
superconducting bound states to the two-particle (Andreev) current. The traditional view of subgap current transport 
in proximityj-NINS and SINS structures considers single-particle tunneling into bound states in the normal region of 
the INS welluEle, implicitly assuming that the normal region of the INS well plays the role of equilibrium reservoir. 
Such a model is appropriate for low-transmission tunnel junctions with low tunneling rate compared to the inelastic 
relaxation rate. However, transparent mesoscopic structures are in a different transport regime where the bound levels 
are well decoupled fmm the superconducting reservoirs, and where injected quasiparticles escape from the INS well via 
Andreev reflect ion .E3o Resonant two-particle current in quantum NINS junctions has been theoretically studied in 
Refs. |ljjl^. In SNS-iunctions, the quasiparticles may undergo multiple Andreev reflections (MAR) before they escape 
into the reservoirs.EE Jpr-a, number of recent experiments with ballistic SNS devices fabricated with high mobility 2D 
electron gas (2DEG)tZrE2l the highly coherent MAR transport regime has been realized. 
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The purpose of this paper is to develop a theory of coherent multiple Andreev reflections which will be applicable 
to 2DEG-based SNS junctions. In 2DEG devices the separation of the superconducting electrodes L is typically larger 
than 200 nm, which is of the same order of magnitude as the superconducting coherence length, £ = ?iw_f/A (vp 
is the Fermi velocity of the 2D electrons), and superconducting bound states are formed well inside the energy gap. 
The presence of bound states in the junctions of finite length gives rise to resonances in the MAR transport, which 
dramatically affects the subgap current. Furthermore, it is possible in 2DEG devices to use electrostatic gates to 
reach the quantum transport regime with a small number of electron modes and_ variable transmissivity. A theory of 
coherent MAR has earlier been developed foi-short superconducting junctions,EIrE3 L £q, where superconducting 
bound states do not play any significant role.E3 Such a theory is consistent with the physical situation in atomic-size 
superconducting point contacts,EZTLHI where quantization of conduction modes has turned out to be very helpful for 
detailed comparison between theory and experiment. The purpose of our study is to investigate the interplay between 
superconducting bound state resonances and coherent MAR in long SNS junctions, L > £o> in the quantum transport 
regime. j-,,-. 

In a number of publications, the coherent MAR approach has been applied to long SNS junctionsSt! 2 ! However, 
these studies were restricted to fully transparent junctions where the bound states are strongly washed out and the 
resonances are not pronounced (in fact, as we will show, at zero temperature the current in such junctions does not 
show any structures). We will study junctions with arbitrary transmissivity, < D < 1 and pay special attention to 
the low-transparency limit, D< 1, where the resonance effects are most pronounced. 

The paper is organized as follows. In Sec. || we derive a ID model for a gated ballistic 2DEG device with one 
transport mode. In Sec. Ill we construct a scheme for calculating MAR amplitudes in terms of wave propagation in 
energy space. In Sec. IV, single current resonances are studied, and Sec. [v| is devoted to a discussion of the interplay 
between resonances in multi-particle currents. The properties of the total subgap current is discussed in Sec. VI. 



II. ID MODEL FOR QUANTUM SNS JUNCTIONS 

We consider an SNS junction similar to the one discussed by Takayanagi et al.0 schematically shown in Fig. |l|. The 
junction consists of a normal conducting channel fabricated with a high mobility 2DEG, which is confined between 
superconducting electrodes. The distance between the electrodes is comparable to the superconducting coherence 
length and small compared to the elastic and inelastic mean free paths and to the normal electron dephasing length. 
The superconductor-2DEG interfaces are highly transmissive, the transmission coefficient typically exceeding a value 
0.75, and the number of conducting modes in the 2DEG channel is controlled by a split gate. 



gate 




FIG. 1. Sketch of the device: a ballistic 2DEG is sandwiched between two superconducting electrodes (S), and an electrostatic 
split gate creates a quantum constriction (dashed line) where only a few conducting modes are open; rare impurities are indicated 
with (x). 



Under these conditions, electrons ballistically move from one electrode to the other while occasionally being scattered 
by rare impurities or junction interfaces. Under a voltage bias applied to the junction, the transport regime corresponds 
to fully coherent multiple Andreev reflections (MAR). To calculate the dc current we will apply the scattering theory 
approachLj lj generalized for superconducting junctions; see Refs. |2^,|3^ and references therein. 

The normal electron propagation through the junction is generally described by the A-channel scattering matrix. 
By assuming the split gate to select only one transport mode, we will characterize the transport through this mode by 
energy-dependent transmission amplitude, d(E), and reflection amplitudes, r(E), r' (E) (the energy E is counted from 
the Fermi energy). The scattering amplitudes satisfying the unitarity relations, dr* + d*r' = 0, \d\ 2 + |r| 2 = 1. The 
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energy dispersion of the scattering amplitudes will introduce the normal electron (Breit-Wigner) and superconducting 
(Andreev) resonances in the scattering problem. The effect r©f narrow Breit-Wigner resonances on coherent MAR was 
earlier studied by Johansson et alBa and Levi Yeyati et alH3 Here we will focus on the effect of Andreev resonances 
and only consider Breit-Wigner resonances which are wide on the scale of the energy gap. This will allow us to neglect 
the energy dispersion of the junction transparency, D = \d\ 2 «const. However, the scattering phases may depend on 
the energy, which yields the Andreev resonances. In the simplest case, this dependence is a linear function within the 
energy interval \E\ ~ A, and we will write it on the form 

d(E)=d e iaE , r(E)=r Q e lbE , (1) 

where a, b are constant. In this case, the scattering properties of the normal channel are similar to those of a ID NIN 
junction. Indeed, the corresponding ID transfer matrix, 

f(F)-( e ~ taE/d ° r*e^- b ) E /d*\ 

m „-i(a-V)E/j AaElj* I ' \ Z ) 



rQe -*(a-b)E/ dQ e iaE/ d * 
can be decomposed into a product of three transfer matrices, 

T(E) = e-^zLiE/Morp^o-jg-icrzLrE/Ata ^ e -icr z k(E)Lirp e -ia,k(E)L r ^ /g\ 

where a z is a Pauli matrix. The first and the last matrices describe ballistic propagation of an electron, with wave 
vector k(E) = ^/2m(Ep + E)/fi w kp + E/A£q, through the right and left N-regions of an effective junction with 
lengths L r = 6A£ /2 and L t = (a-b/2)A£ respectively (from right to left), and the matrix f = e ^^ F L l f^ e za z k F L r 
describes an effective barrier (I). 

Quasiparticle propagation through the, effective ID SNINS junction is described by means of the time-dependent 
Bogoliubov-de Gennes (BdG) equation,E3 

H {x) A(x)e sgn(,)^ t M \ / u(a , t) \ _ / u{X: t) \ 

A(x)e~ s % n ^ eVt / h -H (x) J \ v{x,t) J ~ 0t \ v(x,t) J ' ^> 

where Hq = p 2 /2m — n + U (x) — sgn(x)ey/2 is the normal electron Hamiltonian, U (x) is the impurity potential and V 
is the applied voltage. The superconducting order parameter A(x) is constant within the superconducting electrodes 
and zero within the normal region — Li < x < L r (see Fig. |^). In the further calculations, the impurity potential is 
described by the transfer matrix T in Eq. (^). The spatial distribution of the applied potential along the channel is 
modeled with a step-like function ±eV/2. In fact, the actual spatial distribution of the potential does not play any 
role in this system: it can be included in the transfer matrix in Eq. (^|), leading to an additional energy- independent 
shift in the scattering phases in the matrix T. As we will see later [comment after Eq. (jlT])], the energy-independent 
phases in the T- matrix do not affect the current, and can therefore be excluded. 

I eV/2 



-eV/2 



-Li 







FIG. 2. Spatial distribution of the superconducting order parameter and the electrostatic potential in the junction. The bold 
vertical line indicates the impurity potential. 

The phase difference between the two superconductors follows from the Josephson relation (4> = 2eV/7i) and 
introduces time dependence into the problem. The superconducting electrodes are considered to be equilibrium 
reservoirs where the quasiparticle wave function is a superposition of electron- and hole-like plane waves, 



±ik e x-i{E±u z eV '/2)t /h 



D ±ik h x-i(E±a z eV/2)t/h 



(5) 
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In this equation, the ± signs in the time-dependent factors refer to the left/right electrode, k e,h (E) is the wave vector 
of electron/hole- like quasiparticles, and u(E), v(E) are the Bogoliubov amplitudes. The ratio of the Bogoliubov 
amplitudes equals the amplitude of Andreev reflection for particles incoming from the neighboring normal region, 



l-„(jr\-K E - *9 <E)y/E* ~ A 2 ) /A, \E\ > A 
u 1 ' \ (E-iVA 2 -E 2 ) /A, \E\<A 



(6) 



Since the time dependencies of the wave functions in the two reservoirs are different, the quasiparticle scattering by 
the junction is inelastic, and one has to consider a superposition of plane waves with different energies in order to 
construct scattering states. 

III. CALCULATION OF CURRENT USING SCATTERING STATES 

A. Recursion relations for MAR amplitudes 

We will now proceed with-the construction of recurrences for the scattering amplitudes following the method 
suggested by Johansson et al.Ea. To this end we introduce the wave functions in the left/right normal region (l/r) of 
the junction with respect to the position of the impurity. A particular scattering state, labeled with the energy E of 
the incoming quasiparticle, will consist of a superposition of plane waves with energies E n = E + neV, where n is an 
integer, — oo < n < oo, 



J, 1 p iK 



+ 4+ e_ 







■cb l e-< x 



(7) 



i(E n -<j z eV/2)t/h 



n— — oc 



J,r ik'lx I A,r -ik'lx 



The normal electron/hole wave vector k^ h is here defined as k^ h = k(±E n ), k(E) — yj2m(Ep + E)/fi w kp + E/hvF- 
The meaning of the labels for the scattering (MAR) amplitudes c„ will be explained below. 

Continuity of the scattering state wave function across the left and right NS interfaces determines the relation 
between the electron and hole amplitudes in the vicinity of each interface, 



n ^ 0, a n = a(E n ), 



(8) 



which describes elastic Andreev reflection (l/r indices are omitted). It is convenient to consider scattering amplitudes 
near the impurity (at x = ±0) rather than at the NS interfaces, and to rewrite Eq. (||) for such amplitudes, combining 
the amplitudes of the ballistic propagation through the normal regions with the Andreev reflection amplitude. Then, 
in vector notation, 



Cn± 



(9) 



the modified relation M) takes the form 



U n c n -, n/0, 



(10) 



where 



a n 
a' 1 



(11) 



The phase ip n = 2E n L^ r / A^q — arccos(i? I j/A), characterizing U n , is real inside the energy gap, \E n \ < A, where it 
describes the total energy dependent phase shift due to ballistic propagation and Andreev reflection. Outside the 
gap, \E n \ > A, the phase ip n has an imaginary part which describes leakage into the superconducting reservoirs due 
to incomplete Andreev reflections. 
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By matching harmonics with the same time dependence in Eq. (R), we derive a relation between scattering ampli- 
tudes at the left and the right side of the barrier: 



L (n+1)- — 1 C n+> 



c (n+l)- ~~ 1 c n+) 



(12) 



where the effective barrier transfer matrix T is defined in Eqs. (^), (||). 

The recursion relations in Eqs. (|l^) and (12) couple the scattering amplitudes c n ± into an infinitely large equation 
system. This equation system describing coherent MAR is illustrated by the MAR diagram in Fig. 0. 




FIG. 3. Scattering state in energy space: coefficients c£± correspond to the part of the scattering state at respective energy, 
location and specific direction. Electrons are indicated with full lines and holes with dashed lines; the arrows above (or below) 
each coefficient indicate the direction of the quasiparticle motion. 



The electron part of the quasiparticle injected at the left NS interface propagates upwards along the energy axis, 
the amplitudes for this propagation being labeled with c^. At the injection energy E — Eq (amplitude cj, ), the 
quasiparticle is accelerated across the barrier (I), where the potential drops. Thus, it enters the right normal part 
of the junction with energy E\ (c}_), undergoes Andreev reflection and goes back as a hole (c} + ), entering the left 
normal part of the junction having been accelerated to energy E^ (cj_), and is then again converted into an electron 

(c| i ). The ± indices label the amplitudes after (+) and before (— ) the Andreev reflection for propagation upwards 
along the energy axis. There is a similar trajectory of injected holes, which descends in energy, with the MAR 
amplitudes labeled with c 1 . Due to electron back scattering at the barrier, the upward and downward propagating 

T T I 

waves are mixed, e.g. c 0+ being not only forward scattered into c{_, but also back-scattered into Cq + , which opens 
up the possibility of interference. 

Injection from the left reservoir, shown in Fig. [| generates a MAR path which only connects even side bands at 
the left side of the junction with odd side bands at the right side. Injection from the right reservoir will generate a 
different MAR path, with even side bands at the right side of the junction, i.e. the diagram in Fig. || will effectively be 
mirrored around the barrier (I). Thus, there are two independent equation systems for the MAR amplitudes: injection 
from the left and from right. The l,r labels in the MAR amplitudes can then be omitted since they are uniquely 
defined by the source term and the side band index. 

The transport along the energy axis generated by MAR, from energy E to E ni is conveniently described by the 
effective transfer matrix M n Q, 

c„_ = M n0 c 0+ , n > c n+ = [Mon]" 1 ^-, n < (13) 



M nm = Tn-xUn-i . . . U m+1 T m , n>m, (14) 

where = T^ 1 and T^+i = T for the injection from the left (for the injection from the right, the even and odd side 
band indices are interchanged). For paths within the superconducting gap, \E n ,E rn \ < A, the matrix M„o satisfies 
the standard transfer matrix equation, M nm a z M\ m — <r z , which provides conservation of probability current along 
the energy axis, 

in± = c f n± cr z c n± , f n± =f m±1 \E n>m \ < A. (15) 
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An important consequence of the coherence of MAR is the possibility of transmission resonances in energy space. 
From the form of the M-matrix, 

M n0 = ...f- 1 e lrT ^ k f..., (16) 

it is evident that when tpf, = mir, the two matrices X 1-1 and T will cancel each other, and the probability of transmission 
through this part will be unity, which leads to resonant enhancement of MAR. The solutions of the resonance 

equation, 

(m) ^E^Ll r E]~ 

tp k = fk — ran = — — — arccos — rrnr = 0, (17) 

coincide with the spectrum of the de Gennes-Saint- James levels localized in INS quantum wells!. The corresponding 
bound states are located either on the left or the right side of the junction. 

Without loss of generality, the calculations can be performed with a real matrix T. The transformation to such 
a real matrix is given by T — > V1TV2, with diagonal unitary matrices V± t 2 whose elements are constructed with the 
scattering phases, which are energy-independent. It is clear from Eq. ( |T6| ) that since these energy-independent matrices 
commute with the matrices U n , they cancel each other, and the matrix M„o undergoes a similar transformation. This 
will lead to an overall phase shift of the scattering state which does not affect the current. 

It is interesting to consider the special case of fully transparent junctions, D = 1, which has been studied in the 
literatureEU'LJ. In this case, all matrices T n in Eq. ( |14[ ) are equal to the unity matrix, and the M-matrix takes the 
simple form M„o = exp(icr z Y] —1 ip m ). The length of the junction then enters only through the phase of the MAR 
amplitudes, which drops out of the side band current. Thus the dc cuxijeot of fully transparent SNS junctions is 
independent of length and equal to the current in quantum constrictionsEjtHI. In particular, at zero temperature this 
current does not show any structures in the subgap current. It is also worth mentioning that in this particular case 
of fully transparent SNS junctions, the M-matrix is diagonal and therefore a closed set of recursive relations can 
be derived for the MAR probabilities (not just for the MAR amplitudes, as in the general case), eouivalent to the 
equations for distribution functions derived in the original paper by Klapwijk, Blonder and Tinkhamll3. 

Equation ( |i"3| ) describes "source-free" propagation along the MAR ladder. To complete the set of equations for the 
MAR amplitudes we need to take into account quasi-particle injection, which introduces a source term in Eq. (|l3|). 
To this end, let us consider a quasiparticlc incoming from the left superconducting electrode with energy E, having 
a wave function of the form 



9f(E) = e 



-iEt/h-icr z eVt/2h 



W*" [ U \+ 6 vh e-^ J 



(18) 



The two terms in this equation refer to electron- (v = e) and hole-like {v = h) injected quasiparticles. We now include 
this wave function into the continuity condition at the NS interface at energy E, which gives us the following relation 
between the MAR amplitudes cq + and co_, 

c 0+ = f> c - + Y (19) 



For quasiparticles injected from the right, a similar equation holds with the substitutions e — > h, and Li — > L r . 
Equations ( |l3| ) and (|l^) give a complete set of equations for the MAR amplitudes with the boundary conditions 
£±00 = at infinity. 



B. Calculation of MAR amplitudes 



A formal solution of Eqs. (13) and jl9|), which is useful both for numerical calculations and analytical investigations, 
can be constructed by reducing this infinite set of recursion relations to a finite set by representing the MAR process 



above E n and below Eq by boundary conditions involving reflection amplitudes r n+ and ro— , defined as = c, 



and cj_ = Cq_t _. This gives the following representation for the vectors in Eq. (j|), 



2n+ = 4+ ( r \ + ) , co- = c l _ ( r °- j . (20) 
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The reflection amplitudes r n+ and tq- are independent of the injection, in contrast to the coefficients c n+ , Cq_. 
Furthermore, they are determined by the boundary conditions c±oq = and can be expressed in terms of the matrix 



elements of Mjv„ and M, 



'O(-AT) 



where N 



oo, 



lim M Nn 



N 



lim [M 0( _ N) Y 



ro- 
1 



0. 



(21) 



In other words, the vectors in Eq. (21) are equal to the asymptotical values of the eigenvectors of M-matrices 
corresponding to the eigenvalues which decrease when N goes to infinity. The advantage of introducing the reflection 
amplitudes r n+ and ro_ is that although they have to be calculated numerically, the recurrences which they obey 
do not contain resonances and converge rather quickly. This is in contrast to the matrix M n o, which does possess 
resonances, but which can be calculated analytically in a st raig ht forward way for any given n. 

The solutions of the recursion equations Eqs. (|10|) and (|l2|) can now be explicitly written down. For any given 
energy E we get four different sets of solutions for four scattering states including electron/hole injection from the 

mil m i2 
^ m 2 i m 2 2 

solutions for injection from the left (n > 0) have the form, 



left and the right. Using the formal expression in Eq. (EOT) and the matrix elements of M. 



ii o 



the 



u(l - al)e lip " [S ve + a r -6uh] 



m 2 2 + m 2 ir _e 



2iip 



mi2r n+ e 2liPn — mnr ^r n+ e 2llf ' 0+2llf ' n 



1 



(22) 



The solutions for injection from the right can be found through interchanging e <-> h and calculating all quantities 
with respect to injection from the right. The solutions for n < are calculated in a similar manner. 



C. Calculation of current 



Now turning our attention to the current, we calculate it in the normal region next to the barrier, using the wave 
function in this region, 'J', and assuming quasiparticlc equilibrium within the electrodes. The current then takes the 
form, 

/(,) " ii C dE ( " 2 - ° 2yl £ Im } to ' ,h SB? • (23 » 

e/h,l/r 



where (u 2 — v 2 )^ 1 = \E\/y/E 2 — A 2 = \E\/£ is the superconducting density of states, and the sum is over the four 
scattering states at a given energy E associated with the electron- and hole-like quasiparticles (e/h) injected from the 
left and right (1/r). The current can be divided into parts with different time dependence and expressed as a sum 
over harmonics: 

I(t) =Y,lNe 2meVt/h (24) 

N 

Focusing on the dc (N = 0) component, and calculating the contribution of each scattering state at the injection side 
of the junction, we express the current spectral density J(E) through the probability currents of electrons and holes 
at energies E 2n (Fig. §), 

A 



Idc = r [ dE J(E), 



I E I 

J(E)= ~T~ ( d 2n-^c 2n _ + cl n+ a z c 2n+ J . (25) 

e/h,l/r n=—oc 

These currents coinside with the probability currents j^±, Eq. (|l5|) , flowing along the energy axis. 

It is convenient to introduce a leakage current J n , defined as the difference of the probability currents before and 
after Andreev reflection, 

e/hd/r ? 



7 



J„ represents the amount of probability current from all the scattering states injected at energy E and leaking out 
of the junction at energy E n (Fig. |^). The leakage current is zero inside the energy gap due to complete Andreev 
reflection, J n = 0, \E n \ < A [cf. Eq. 

The explicit expression for the leakage current for n^O follows from Eq. (^6|) after insertion of Eqs. ( ^2[) and (|Tc|), 

Jn = y (1 - |a | 2 )(l - K| 2 )(l + |ro-a | 2 )(l + |r„ + a»| 2 ) 

77^ 1^22 + e 2lv3 °ro-TO2i — e 2lv,i r„ + TOi2 — e^oe 2 *^"^— fn+Tniil 

1 1 r 

It follows from Eq. (|2^) that the leakage currents are positive for all n ^ 0, J„ > 0. One can also show that they 
satisfy the inequality X)n^o ^ n — ^' wn i cn is a consequence of the conservation of probability current: the leakage 
current of all side bands except of the side band n = does not exceed the probability current injected into four 
scattering states. Furthermore, the leakage current satisfies the important detailed balance equations, 

J_ n (£) = J„ (£_„) (28) 

i.e. the leakage at energy E_ n due to the injection at energy E is the same as the leakage at energy E due to injection 
at energy E- n . 

Using the continuity of current across the barrier, = j^ n+1 ^_, guaranteed by the transfer matrix T, we can 
express the probability currents in Eq. ( |25| ) through the leakage current, Eq. (|26|), 

I F\ °° 

e v^-=e j ^ ™>° ( 29 ) 

e/h,l/r k=n 



E ^'• ? ™+ = ~ E J - fc ' n < °' 

e/h,l/r k— — n 

by adding and subtracting consecutive terms in the sum. The spectral density of the dc charge current Eq. ( |25| ) can 
then be written on the form 

J(E) = J2nJn(E), (30) 

n 

since J n appears in n probability currents. This formula has a clear physical meaning: the contribution to the charge 
current of the n-th side band is proportional to the leakage current of the side band times the effective transferred 
charge ne. 

The detailed balance of the leakage currents, Eq. (p8|), allows us explicitly to prove that at zero temperature the 
scattering processes between (occupied) states with negative energies, E, E n < —A do not contribute to the current, 
in agreement with the Pauli exclusion principle. Indeed, by separating the contributions from side bands with n < 
and remembering that the leakage current is zero within the gap, we get for zero temperature, 

dEJ2nJ n (E) = V — / dEJ n (E) + / dEJ n {E) - / dEJ_ n (E) , (31) 

n^O n>0 11 V-™ J-oo J 

where the first and the third terms cancel each other by virtue of Eq. (|2|). At finite temperature, these two terms 
produce current of thermal excitations while the second term gives the current of real excitations created by the 
voltage source. Keeping only this term which dominates at low temperature, we finally get, 

r-A 

Ida = E J »' 7 " = ^{neV - 2A) / dEJ n {E) tanh(|£|/2£; B T). (32) 

„ >0 JA~neV 

We end this section by noting a technically useful symmetry in the current density, namely J„(E) = J n {—E — neV), 
seen from the explicit form of the M n o-matrix. This allows us to reduce the integration interval in Eq. ([32]) to 
-neV/2 < E < -A. 
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IV. CURRENT IN TERMS OF n-PARTICLE PROCESSES 



The approach formulated above provides necessary foundations for numerical calculation of the current for arbitrary 
transparency and length. However, to get a full understanding of the rich subgap structure in the current-voltage 
characteristics, which may seem quite random, especially for intermediate transparencies and lengths (see Figs. |l4]-[l6|), 
we will conduct a detailed analytical study of the limit of low transparency D < 1, The separation of currents into 
n-particle currents, Eq. (j32|), is our basis for analysis, and we will study each current /„ separately. 

As explained in the previous section, the de Gennes-Saint- James levels, Eqs. ( |l6| ) and (|l7|), are important for the 
current transport through the junction, leading to resonant enhancement of the current. Our main attention in this 
and the next section is on the calculation of the position, height and width of the main current peaks and oscillations 
which have the magnitude of order D. To simplify notations, left/right injection indices are omitted in most cases. 



A. Single-particle current 

The single-particle current, which dominates at large applied voltages, has, according to Eq. ([32]), an onset at 
eV = 2A. The full numerical solution for the single-particle current is plotted in Fig. |[ The current shows pronounced 
oscillations, and the magnitude of the slope at the current onset strongly depends on the junction length. 
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FIG. 4. Single particle current for symmetric junctions Li — L r — L/2 for different junction lengths; the junction trans- 
parency is D = 0.1. The current onset for the short junction (L = 0) disappears for junctions with finite length (bold line); for 
L = nn^o, the onset appears being roughly n + 1 times smaller than the onset for L — 0. 



To understand this behavior, we analyze Eq. ( P7| ) in the limit of small transparency D < 1, i.e. in the tunnel limit. 
First we note (see appendix |A|) that the reflection amplitudes r n+ and ro- may be expanded as 

r n+ = (~l) n VR + 0(a 2 n+1 D) (33) 



r _ = y/R + 0(a 2 _ 1 D), 

After inserting the explicit form of M%o = T together with the expansion ( |33| ) into Eq. (|27|) and putting R = 1 , we 
can write the single particle current on the form 

en r~ A 

h = ~rO{eV - 2A) / dE [N l {E)N r {E 1 ) + N r {E)N l (E 1 )] t&nh(\E\/2k B T), (34) 

" JA-eV 

where 



N"(E) = \EWE^ . (35 ) 

E 2 - A 2 + A 2 sin 2 {2EL Ur / A£ ) 



In analogy with the tunnel formula for the current,! N l ' r is identified as the tunneling density of states (DOS) on the 
left/right side of the junction. In Fig. ||the energy dependence of the DOS is presented. The deviation of this DOS 







from the normal metal density of states is a manifestation of the proximity effect. The expression ( |35| ) for the DOS 
has earlier been derived for proximity NS sandwichesEJo'tj. Note that the DOS in our case is constant throughout the 
N-regions. In junctions with arbitrary length, the DOS usually approaches zero at the gap edge \E\ = A (Ref. [To"| ). 
Exceptions are junctions with lengths L^ r = m7r£o/2 where a bound state splits off from the gap edge. In this case, 
the DOS diverges at the gap edge. The quantum well structure of the SNS junctions also give rise to quasi-bound 
states in the continuum spectrum, \E\ > A, seen as oscillations in the DOS. 




FIG. 5. Density of states in the N region for different lengths L of the region. Singular behavior of the DOS for short 
junctions, L = (equal to the DOS in a superconductor), is suppressed for finite length junctions. The amplitude of the 
first oscillation increases as the length increases, indicating accumulation of the spectral weight at the energy gap edge and 
formation of a bound state for L = 7r£o- 

The single particle current in Eq. |34|) is written as the integral over the product of the DOS at the entrance 
energy E and the exit energy E + eV. The latter depends on the applied voltage, as well as the integration interval, 
and therefore the DOS oscillations produce oscillations of the current I\ as a function of voltage (Rowell-McMillan 
oscillationsotj). The oscillations become more pronounced when the junction is sufficiently long, and the differential 
conductance may even become negative. It is also clear that the DOS oscillates as a function of the length of the 
junction, which give rise to oscillations also in I\. 

In short junctions, L <C £o, the current onset at eV — 2 A is very steep, see Fig. [|. In junctions with finite length, 
the current onset is smeared and replaced with a smooth oscillating behavior. This can be directly related to the 
smearing of the singularity in the DOS at the gap edge. The length where the crossover between these two behaviors 
occurs can be taken as a measure of when finite length effects become important. To estimate this length, we write 
equation (|27|) for small lengths L <C £o, near the threshold, eV = 2 A + ft, ft <C A, keeping the first order terms in D 
in the denominator. For a symmetric junction, Li = L r = L/2, we get 

eAtanh(A/fc B T) f w D sin 2 z dz 

h = = ' " daV da^ (36) 



4ft / v 4ft 

From this formula it is clear that for short junctions (L = 0), the current onset has the width ft ~ AD/4. If L is 
of the order of £oV^D/2, the size of the onset has substantially diminished, and there is no visible onset at eV = 2A 
when L 3> £,q\/~D/2. This crossover between steep onset and smooth behavior, which happens already for quite short 
lengths if D is small, can be interpreted in terms of a bound state which is situated exactly at the gap edge in short 
junctions (L = 0), and which moves down into the gap when L > 0, the effect becoming fully pronounced when the 
distance from the gap edge, Hvf/L, exceeds the dispersion of the Andreev state, \/~DA 1 in symmetric junctionsEl. 

When L^ r approaches 7r£o/2, the lowest quasi-bound state in the continuum spectrum approaches the gap edge. 
This leads to an accumulation of the spectral weight at the gap edge and reappearance of the singularity in the DOS, 
which results in the reappearance of a sharp current onset at eV = 2A, but with smaller magnitude; see Fig. ^ 
(Li t r = 7t£ )- 

It is of interest to note that in our calculations, based on the scattering theory approach, the bound states are 
not directly involved in the single-particle transport, which therefore is non- resonant and shows no subgap resonance 
peaks. Within the tunnel model approach the situation is qualitatively different: the DOS in Eq. usually includes 
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the contribution of the broadened bound states, and therefore the single-particle current exists and has pronounced 
resonant features at subgap voltages eV < 2A. This difference results from the fact that, within the tunnel model 
approach, the superconducting bound states are implicitly assumed to be connected to the reservoirs (broadening 
due to inelastic interaction), which allows a stationary current to flow through the bound states. In contrast, within 
the scattering approach, the bound states are disconnected from the reservoirs and have zero intrinsic width. In this 
case the bound states obtain their width only due to higher order tunneling processes involving Andreev reflections, 
which are manifested by the resonant multi-particle currents. In practice, the relevance of the multi-particle versus 
single-particle mechanism of the subgap curpnt transport is determined by physics and depends on the ratio of the 
corresponding dwelling and relaxation timesc3. In this paper, the inelastic relaxation time t% which determines the 
width of the single-particle resonances is assumed to be much larger than the dwelling time of the most important 
two-particle current, Tj ^> hvp/LD. 



B. Two-particle current 

The two-particle current I2 in quantum point contacts, (L -C £0) is of order D 2 when eV < 2A and of order 
D 2 h\D when eV > 2 A (Ref. |25|). For finite length junctions, the situation is different. For the MAR paths where 
the energy of the Andreev reflection coincides with a bound state, the current spectral density j\ is of order unity, 
due to resonant transmission through this state. For low transparency D <C 1, this gives a sharp concentration of 
the current density around the resonant energies. In this limit, the two-particle current is well described by the sum 
of contributions from these resonances, and to evaluate them we examine the energy dependence of J 2 close to the 
resonant energies, E\ = + SE. Let us consider the contribution to the leakage current [J2] 1 from quasi-particles 
injected from the left. As shown in appendix [TJ, in this case Eq. ( |27| ) reduces to the standard Breit-Wigner resonance 
form 

r (m) r (~^ 



= i V^r— ; tt~2 ' (37) 




where the tunneling rates T n are given by li m) = N l {E n )D/2^ m \ n = 0, 2, and 

2L r A 



' dE 

and the position of the resonance is shifted by 



(38) 



SE ™ = -^Im \ i^— + \ . (39) 



An analogous result is valid for quasiparticles injected from the right. 

After integrating over energy, the two-particle current in the resonance approximation may be written on the form 

HeV) =YY^6(eV- - A) 2 * DA ^(E^-eV)N*(E^+eV) 



i—l.r m>0 



where the summation is over the positive bound level energies, < E^ < A, and the DOS iV* should be calculated at 
the injection side of the junction and / (m) (T, V) = (1/2) [tanh ({eV - £:( m ))/2fc B T) + tanh {{eV + £( m ))/2fc B T)] . 
According to Eq. ([i7j|), the two-particle current /2(eV) increases in a step- like manner in the voltage region A < 
eV < 2A. The steps occur at every voltage where a new resonant channel through a bound state opens up, at 
e y{m) = A + E( m \ We note that the step positions depend on temperature and approximately scale with A(T). 
Each current step has the height of order D. As seen from Eq. ©, the contribution to the current of a particular 
bound state E^ is modulated, as a function of voltage, by the oscillations of the density of states at the entrance 
and exit energies, N(E^ m ' ) ± eV). In other words, the pronounced oscillations of the two-particle current seen in Fig. || 
reflect how close the entrance and exit energies E^ ±eV are to a quasi-bound state in the continuum. For eV > 2 A, 
the two-particle current I2 oscillates around a constant value with an amplitude of oscillation decreasing as A 2 /eV 2 
for large voltages. 

It is interesting j-tp compare the resonant structures of the two-particle current with the resonant structures in 
NINS j unctions .Ej'Ej In NINS junctions, the resonant current steps occur at eV — E^- m \ and they do not have any 
modulation because the DOS on the normal side of the junction is constant. 
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FIG. 6. Two-particle current in symmetric junctions Li — L r — L/2 for different lengths; the junction transparency is 
D = 0.1. The resonant process shown in the inset becomes possible when eV > A + |_E^ m - ) |. 

The distance between the resonances and the resonance widths are proportional to the bound level spacing, and 
they decrease in long junctions. For sufficiently long junctions, the two-particle current may thus give the appearance 
to include a series of peaks, as shown on Fig. [?]. In symmetric junctions, the bound state energies at both sides of the 
barrier will coincide, reducing the number of steps with a factor of two, and giving current steps of double height. 

We will conclude this subsection by noting that the difference between the full numerical calculation of the two- 
particle current and the resonant approximation given in Eq. (|40| ) is rather small already when D = 0.1, as can be 
seen in Fig. ^. 



0.5 




FIG. 7. Comparison between the approximate expression for the two-particle-current in Eq. (^0|) (solid curve) and the full 
numerical solution (crosses): D — 0.1, Li — L r = L/2 = 5£q- 

C. Excess current 

Excess current in SNS junctions, i.e. the difference between the current in superconducting junction and in the 
normal junction at large voltage, 

I exc = I-G N V + 0(A/eV), (41) 

is commonly considered as a measure of the intensity of Andreev reflection. In tunnel SIS junctions and low- 
transmissive point contacts the excess current is siHaJl, I exc « D 2 eA/irh, D<1, while in fully transparent contacts 
the excess current is large, I exc = 8eA/37r?i, D = l.Ea Accordingly, one would expect large excess current in long SNS 
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junctions due to the resonant enhancement of the two-particle current. However, the excess current is small because 
of a large deficiency, of order D, of the single-particle current caused by the broadening of the current onset at the 
threshold. As we will show, the single-particle and two-particle currents undergo a fine cancellation, yielding small 
net excess current of order D 2 when D <C 1 . 

The excess current has contributions only from the single- and two-particle currents, since all higher order currents 
include at least one Andreev reflection outside the gap whose probability is of order (A/eV) 2 . In the limit of large 
voltage, eV A, the relevant part of the current in Eq. (p2h then takes the form 



h 



4De 



-A 



dE 



iV/2 



(l-a 2 )(l + Ra 2 ) 
R 2 a A - 2R Rc{e 21 ^} 



(42) 



8D 2 e 



-A 



dE 



eV 



«*1 



l + i? 2 |ai| 4 -2i?Re{e 2 ^i}' 



These equations are written for symmetric junctions, L r = Li = L/2, and for zero temperature; small Andreev 
reflection amplitudes < |a(eV/2)| <C 1 have been neglected in Eq. (p7[). The behavior of the current in Eq. (^2|) as a 
function of voltage is presented in Fig. || for different lengths. It is clearly seen that the limiting value of the excess 
current is approached much faster in finite length SNS junctions compared to point contacts (L = 0). In Fig. || the 
excess current behavior with respect to the junction length is presented for different transparencies. 



0.6 
0.5 

< 

I 0.4 



0.3 

+ 
O 



0.1 





.... L=0 




— L = 2 ^o - 




... L=7^ Q 


■ 
i 
i 
l 

i 

* 

\ 
\ 

\ 

\ 

X 





eV/A 



8 



10 



FIG. 8. Deviation of the current from its asymptotical value at V = oo: the excess current value is approached much faster 
in finite length junctions, here shown for D = 0.3. 
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FIG. 9. Dependence of the (normalized) excess current on the junction length for different transparencies. For fully transpar- 
ent junction, D = 1, the excess currents are identical for all junction lengths; the excess current increases for small-transparency 
junctions. 
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To analytically examine the excess current in the limit of small transparency D <C 1, it is convenient to start with 
equations (34) and (^p|). To first order of D the excess current assumes the form (T = 0), 



— 1 l 



jexc 
1 2 7 



(43) 



jexc 
1 1 



AeDA 2eD 



/ [N l (E) +N r (E) - 2] dE, 

J A 



2 2— 1 hn( m ) 

l/r, m>0 



Let us consider the contributions to the single-particle current from the left electrode, 

\iri = - 2 -^ + ^ l°° i N \ E ) - 1] dE - ( 44 ) 

n h J A 



Inserting N (E) from Eq. (p5[), this equation can be transformed to the form, 

, = 2_eD [°° f Et E\ dE = e_D [™ jin^^ 

h Ja \e + ^ 2 sm 2 (2EL l /A^ ) £j h J ^ V + S m 2 (2ELi/A£o) ' 



where ^ = ^/E 2 — A 2 . It is now possible to analytically continue the integral in the upper half plane which will reduce 
the integral to a sum over the residues of the poles given by the equation £ 2 + sin 2 (2i?i;/A^o)- Comparing this 
equation with Eq. ( p7| ) we find that the poles coincide with the energies of the bound states in the gap. The excess 
current contribution from the left-injected single particle current is thus 

Kl-Tsi-rf, <«> 



m>0 



where \l^ c \ is the contribution to the two-particle current from the bound state resonances at the left electrode. 
A similar relation is derived for current from the right electrode. Thus, there is exact cancellation of the excess 
single-particle and two-particle currents to first order in D. 

It is interesting to note that the cancellation effect is related to the conservation of the number of states in a 
proximity normal metal compared to the conventional normal metal. It follows from Eq. ( p} ) that I\hj2eDA is equal 
to the difference between the number of continuum states in the proximity metal and the total number of states in a 
conventional metal, while, on the other hand, the number of the bound states, is equal to 



/•A pA 

dE ${<P{E) -mir)= dE V S(E - eK™))/?/™) = I 2 h/2eDA, (47) 



m>0 u m>0 



according to Eq. (|43|). 



V. INTERPLAY BETWEEN RESONANCES 



For processes with several Andreev reflections (n > 3), the possibilities for resonances increase. Every Andreev 
reflection energy may coincide with a bound state energy and thus be resonant. For some specific voltages, more 
than one resonance is important, creating a situation of overlapping resonances, which can enhance the current giving 
peaks in the current-voltage characteristics at these voltages. 



A. Three-particle current 

The three-particle current I3 has a non-resonant value of order D 3 . However, ^3 is enhanced to order D 2 when the 
energy of one of the two Andreev reflections coincides with a bound state energy. 
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FIG. 10. Three-particle current in symmetric junctions Li — L r — L/2 for different lengths; the junction transparency is 
D = 0.1. The MAR path with two overlapping resonances, shown on the inset, generates a current peak with height proportional 
to D. 



For the applied voltage equal to the difference between two bound state energies, eV^ km ^ — E^ — E^ k \ two 
resonances occur simultaneously, i.e. form a resonance consisting of two overlapping single resonances; see the inset 
in Fig. [l(]. This will enhance the current to order D close to this voltage, giving a peak in the IVC. The number of 
peaks is equal to the number of bound state pairs. The peaks are located in the voltage interval 2A/3 < eV < 2 A; 
we note that the peak positions are weakly dependent on temperature. 

To evaluate the height and the width of these peaks, we study the contribution from overlapping resonances at 
Ei w E^ < and at E 2 « £ (m) > 0. Close to these energies, the phases ip[ k ^ and , defined in Eq. (|l7|), are 
close to zero, and we find the current spectral density for injection of a quasiparticle from the left (see Appendix^), 



[Jg](*"0.' 



D 3 N l (E)N r (E 3 ) 



D 



^ {k \ 2 m) + iD ^[ k) Nr(E 3 ) + ^N\E) 



(48) 



We now expand ip[ k \ </4"^ 
voltage, SV = V - V {km ^ 



in the deviation from perfect overlap in energy, 5E = (E\ — E^ + E 2 — E^ m ^)/2, and in 



and find, using D <C 1, from Eq. (48) 



[ME)} 



(km),l 



(fe)-p(m) 



dtx ' r 



(5E+SE-/A 2 - D/^WjjM) + A 2 



where A = (T 3 mj 6E + + T {k) SE_)/A, 5E± = SE ± eSV/2. The energy dependence of the current in Eq. ([19) has 
the form of two resonant peaks with width ~ TA ~ DA/r] split by the energy interval ~ s/DA/r) at SV = 0, the 
peak splitting increasing with increasing SV. After integration over energy, the overlapping resonances give a current 
contribution in the form of a current peak (fc^T <gC A), 



(49) 



I. 



(km) 



(SV) 



3De 



7rA 



2N l N 



1 + ?j( fe )ri( m ) 



SV \ 2 r/W N r + r]( m ) N l ' 



DA 



(50) 



In this equation, the densities of states N > r are taken at the entrance and exit energies, N l (2E^-E im) ^j and 

N r (2E^ m ^ — E^), and the temperature is taken to be zero. A factor of 2 has been included in Eq. ( [SO] ) to take into 
account the similar resonant process for injection from the right, where E\ = — ij( m ) and E 2 — ~E^ k \ 

The curve for the three-particle current versus voltage thus consists of peaks with heights of order D and half- width 
Ty = \/DA/rj on top of a background of order D 2 . The background current increases with voltage in the interval 
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2 A/3 < eV < A as more single resonances come into the integration region. In the interval A < eV, the background 
current decreases due to broadening of the resonances because of leakage associated with incomplete Andreev reflection 
outside the gap. 

In long symmetric junctions the current peaks form an interesting triangular pattern. To see this, we first note that 
if the bound state spectrum were perfectly linear, several of the peaks described by Eq. (|48|)- (l50| ) will be situated at 
the same voltage since eV^ km '> = E^ - E^ = £'( m + 1 ) - (see also the inset in Fig. [□]), and thus the total 

number of peaks will be reduced while their respective height will be increased. Since the bound state spectrum is not 
linear, the peaks show splitting. However, the deviation from linearity is small and in practice the peaks form clusters, 
giving combined peaks with height roughly equal to the number of clustering peaks. In the interval 2A/3 < eV < A, 
the number of peaks in a cluster increases in steps of three from 1 to 4, etc., up to the number of bound states. In 
the interval A < eV < 2A the number of peaks in a cluster decreases in steps of one. This gives an appearance of a 
"peak triangle" for very long junctions, shown in Fig. [n]. 




FIG. 11. "Peak triangle" of three-particle current for long junction: Li — L r — L = 10£o, D = 0.1. Every peak of the triangle 
consists of a number of tightly positioned resonances due to nearly equidistant bound state spectrum (resonance cluster). The 
number of resonances in a cluster is, from left to right, 1, 4, 7, 6, 5, 4, 3, 2, 1. The inset shows an example of resonant MAR 
paths forming a cluster. 

This "peak triangle" is further enhanced by the background current, which has a similar triangular form, as 
explained above. 



B. Four- particle current 

The four-particle current has a non- resonant value of order D , which is enhanced to order D 3 when the energy 
of one of the three Andreev reflections coincides with a bound state energy. Similar to the three-particle current, 
overlapping resonances can enhance the magnitude of the current J4 to the order D for those voltages where both 
the first and the third Andreev reflections coincide with the bound states, as shown in the upper inset in Fig. |l2] . 
Indeed, it is clear from the explicit form of Af 40 = f e ia ^ 3 f ~ x "^f e ia ^f ~ x that when ip 1 = kir and cp 3 = rrnv, 
then M40 = (— l) fe + m e l<JlV2 , i.e. the transparency of the MAR trajectory is enhanced to unity. Other combinations 
of the resonances, e.g. when the first and the second Andreev reflection occur at bound state energies, will produce 
peaks of order D 2 or smaller, as described in Appendix 

Focusing on the double resonances which produce large (~ D) current peaks, we find that in short junctions with 
just one pair of bound states, ±E^-°\ the double resonance will occur at voltage eV — E^°\ provided the energy of 
the bound state is within the interval A/2 < E^ < A. The spectraLdensity of the current has a form similar to the 
one in Eq. (|49|), the major difference being the small peak splitting^, ~ DA/ '77. The height of the resulting current 
peak {k B T <€. A) is 

ih]max - h l +[a( 2£(o))] 4 ' (51) 
where a(2E^) is the Andreev reflection amplitude at energy 2E^ \ 
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For longer junctions, there are many possibilities to have overlapping resonances. Two bound states at one side of 
the junction with energies E {k) < and E {m) > can give a peak in 7 4 if (E (m) - E [k) )/2 = eV > A/2. Although 
the height of all peaks is roughly proportional to D, numerically the heights (and widths) of the peaks may vary 
considerably depending on the position of the second Andreev reflection. If the second Andreev reflection does not 
occur at the energy of a bound state, the situation is similar to the one described above; see lower inset in Fig. [T^ . 
However, if a bound state is close to the energy of the second Andreev reflection, then the current spectral density 
IfiE) consists of the three full-transmission peaks with widths ~ DA/77 which are split within the interval ~ s/DA/rj 
(triple resonance). The triple resonance has larger spectral weight compared to the double resonance, which results 
in the larger height and width of the current peak. 

Rigorously speaking, a triple resonance can only occur in asymmetric junctions because it requires equal distance 
between neighboring resonances, while the bound state spectrum in symmetric junctions is not equidistant. However, 
in long junctions, the deviation from the equidistant spectrum is small, and quasi-triple resonances may therefore 
occur also in long symmetric junctions. 




FIG. 12. Four-particle current in symmetric junctions Li = L r = L/2 for different lengths; the junction transparency is 
D — 0.1. The four-particle current in short junctions is not visible on the scale in the figure. The solid-line peak and the small 
dashed-line peaks are due to double resonances, illustrated by the MAR diagram in the upper inset. Large dashed-line peak is 
due to a quasi-triple resonance in the MAR path. An effective four-barrier structure equivalent to this MAR path is shown in 
the lower inset. 

This effect can be observed in Fig. [l^, where the four-particle current for a symmetric junction with length L = 
7£o > 2-7r£o consists of three peaks with different heights: the central peak corresponding to the quasi-triple resonance 
while the two side peaks corresponding to the double resonances with the heights given by Eq. (|5l|). 

Finally, it is worth noting that, similar to the situation for the three-particle current, the peaks will form clusters, 
giving a smaller number of current peaks than the number of pairs of bound states in long junctions. 



C. High order currents 



The studied properties of multiple resonances in three- and four-particle currents allow us to make some general 
conclusions about resonant behavior of the high order multi-particle currents which determine the total current at small 
voltage. The non-resonant magnitude of an n-particle current is of order D n at the threshold voltage, eV n = 2A/n, and 
therefore the total non-resonant current expopentially decreases with the applied voltage (in transparent junctions, 
D ~ 1, the current is exponentially small ato eV < A^/l — D). However, multiple resonances may enhance the 
magnitude of the current by several orders of D. The major question of interest here concerns the maximum value of 
the resonant current, in particular whether it can be of order D at arbitrary small voltage. 

To obtain such large current at small voltage, it is necessary to achieve a transmission probability through a high 
order MAR path equal to unity, which implies that the energy of at least every other Andreev reflection must coincide 
with a bound state (cf. the discussion in the previous subsection). For n > 4, this means that three or more bound 
states must be approximately equidistant in energy. Since the bound state spectrum is non-equidistant, Eq. (|T^), 
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this is generally not possible if the resonances are narrow; therefore, in junctions with arbitrary geometry and small 
transmissivity there are no large current peaks below the voltage eV — A/2. 

However, the possibility of a large resonant current exists for junctions with sufficiently large transparency. To 
find the relevant transparency, let us consider a very long symmetric junction and assume for the moment that the 
bound state spectrum is equidistant, E {m+ ^ - E<- m > = const. Then, from mapping of the n-th order MAR process 
on a ID multi-barrier structure (see Fig.[l3|), it is clear that if the applied voltage is commensurate with the level 
spacing, e.g., eV — £K m+1 ) — E^ m \ the multi-barrier structure is periodic, and full transmission is achieved leading to 
a current peak. This conclusion is valid also for a non-equidistant spectrum if the variation of the interlevel distance 
does not exceed the width of the full-transmission band. The deviation of the bound state spectrum from the best 
linear fit does not exceed the value 0.33A£o/-k, FigjH- On the other hand, the width of the full-transmission energy 
band is ~ y/DA/rj for equidistant spectrum and for n — > oo. Thus one should expect large current structures in 
long symmetric junctions with transparency D > 0.1 to occur at voltages eV > A^q/L. In junctions with smaller 
transparency, large current structures may appear only at eV > A/2, as explained before; see Fig. [l6| It is also easy 
to see that in asymmetric junctions, where the width of the full transmission band for an equidistant spectrum is 
~ DA/n (since the relevant resonances at one side of the junction are weakly coupled to each other through the MAR 
process), large resonant current at small voltage may exist if D > 0.33. Our numerical investigations confirm that in 
symmetric junctions when D is of the order 10 -3 , the multiple resonances are completely blocked and current peaks 
are exponentially suppressed at eV < A/2. 




FIG. 13. Mapping of a high-order MAR path on a multibarrier structure: for an equidistant spectrum, full alignment of 
positions of bound levels (indicated by bold lines) is possible for voltage eV = E^ m+1 ' — E^ m \ yielding a full-transmission 
band. The deviation of the real bound level spectrum from a best linear fit is shown by the thin line. 




FIG. 14. Total current in symmetric junctions Li — L r = L/2, for different lengths; the junction transparency is D = 0.1. 
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FIG. 15. Total current in asymmetric junctions Li = 0, L r = L, for different lengths; the transparency is D = 0.1. 




FIG. 16. Total current for different junction transparencies D. The junction is symmetric, with length Li — L r = 2£o- 



VI. SUMMARY 



Adding up the contributions to the current calculated in this paper, we arrive at a rather complex form of current- 
voltage characteristics (IVC) at subgap voltages, as shown in Figs. |l4[jlq . Nevertheless, the analysis of the tunnel 
limit allows us to classify various subgap current structures. Here we will summarize the results of this-classification. 
As a reference system we will take a short (L = 0) junction where the form of the IVC is well studiedEi The current 
structures in short-junctions can be interpreted as resonant features due to quasi-bound states situated at the edges 
of the energy gap23, the resonant conditions selecting voltages equal to the gap subharmonics, eV — 2A/n. This 
subharmonic gap structure of the short junction gradually changes with increasing junction length as bound states 
move down into the gap, giving rise to IVC structures with steps, oscillations and peaks. Major points are: 

(i) The current in the subgap region is considerably enhanced, compared to the short junction case. This effect is 
present as soon as the effective length L/£o is comparable to, or larger than, the square root of transparency of the 
junction, L/£o ~ y/~D. 

(ii) The main onset of the current in short junctions at eV — 2A shifts downwards in voltage to the value eV = 
A + EW where is the energy of the bound state. This shift is caused by the resonant two-particle current giving 
a contribution to the total current of the order of the single-particle current. 

(iii) For longer junctions, the current onset transforms into a staircase within the voltage interval A < eV < 2 A with 
the number of steps corresponding to the number of bound states, the step positions being given by eV = A + E^l 
This is due to the resonances in the two-particle current transported through bound states. Resonant channels 
open up, one by one, as the voltage increases and bound states enter the "energy window" available for two-particle 
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processes. The current plateaus are not flat but modulated because of oscillations of the density of continuum states. 
The period of the modulation is roughly equal to the interlevel distance, and it decreases with the junction length. 
The amplitude of the modulation, on the other hand, increases with the junction length. Thus, in long junctions, 
the current structures take the form of a series of peaks (see Fig. |l4|) within the voltage interval A < eV < 2A. The 
position of the peaks has pronounced temperature dependence, scaling with the temperature dependence of the order 
parameter, while the distance between peaks has a weak temperature dependence. 

(iv) There is another series of the current peaks whose positions only weakly depend on temperature and are entirely 
determined by the bound state spectrum: eV = — E^ and eV = — E^)/2. These peaks are caused by 
the overlap of two resonances in the three- and four-particle currents and they exist in the intervals of applied voltage 
2 A/3 < eV < 2 A and A/2 < eV < A respectively. The heights of these peaks are comparable with the heights of 
the two-particle current structures (~ D). 

(v) At voltage smaller than eV — A/2 the resonant current structures generally become smaller in magnitude 
(at least by one order in D) if the junction transparency is sufficiently small (D <C 0.1), and the current decays 
exponentially when eV approaches zero (although for some particular junction lengths there could be huge (~ D) 
current peaks caused by multiple resonances). This qualitative difference of the IVC below and above eV = A/2 
allows one to expect a cross over from power to exponential dependence of IVC in multichannel junctions. 

(vi) In transparent junctions, all current structures will persist but become smooth; appreciable current will appear 
below eV = A/2 as soon as D > 1/3. The current structures completely disappear in fully transparent junctions, 
D = 1, where the IVC does not depend on the junction length; see Fig. [16]. _ 

(vii) At voltage larger than 2A, the current undergoes oscillations, similar to Rowell-McMillan oscillationscl, and 
the excess current is approached much faster than in short junctions. In low transparency junctions the excess current 
is small, I exc ~ D 2 ,D <C 1, despite strong Andreev reflection and large pair current 1% ~ D. 
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APPENDIX A: APPROXIMATION FOR r - AND r n+ 

In this appendix, the expansion is derived for the reflection amplitudes in Eq. ( |33|) for a quasiparticle injected from 
the left. From the definition of tq- and r(_i)_, Eq. (p0|), we know 

6o-=c£_( r °f) (Al) 
I ( r (-i)- 



*-u- = i-D- { T" ) • (A2) 

They are related as 

c - = Mq_iZ7_iC(_ 1 )_, (A3) 
where M0-1 = T an d U—\ = e 1(7zip - 1 . From this relation, we find ro- in terms of f*(_x)_ as 

1 + VRr { -i)-e 2i <P-i \ l-x y v ; 

where x = (1 — \ZR)r^_ 1 -j_e 2lip - 1 /(l + r(_ 1 )_e 24V - 1 ). When \x\ <C 1, we can make an expansion in this parameter to 
get to the form 

r - = VR + D 1 }- l) - e ^ = VR + 0{a\D), (A5) 

Similarly we also get 

r n+ = (-1T^+D ^^^l = (-irVR + 0(a? l+1 D). (A6) 
1 + (-l)"r ( „ +1)+ e-^"+i 
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APPENDIX B: RESONANCE IN TWO-PARTICLE CURRENT 



In this appendix, we derive the resonant form of the two-particle current, Eq. <p7\), for a quasiparticle injected 
from the left. The definition of M20 is M20 — Te (TzVl T~' 1 , which using the pseudo-unitarity of the transfer matrices 
<j,T'<j, = T^ 1 can be written on the form 



M 20 = -S= sin <PiT<r, + e- ia ' Vl . 



It simplifies in the limit D<1, |Vi I = Wi — m7T \ *C 1 to 



M 2a = ^-^{2i^ n) {l + a x )+D 



(Bl) 



(B2) 



Inserting the simplified expansion of M 2 q an d the expansion of r 2+ and ro_ from Eq. ( |33| ) into Eq. (27), as well as 
putting R = 1, the leakage current density takes the form 



[Mm 



l-|ao| 4 l-|a 2 | 4 



D 



2i<p\ 



(m) 



D ( 1 + e 2i,p ° 1 + e 2lV2 



2 V 1 — e 2 ^ 1 — e 2i¥>2 
We make an expansion of the phase tp^ = rf m ^{E — E^)/A — rf m ^5Ej A, where 



V 



9 1/3 



2L r 



(») Co ^A 2 - (£(™)) 2 



The two-particle current density now takes a Brcit-Wigner form 



[ME)} 1 



F (m) r (m) 
1 1 2 



-i("0 , p(m) ' 
• > 1 2 



A / V 2 

where the tunneling rates are given by = N l (E , 2 ) D/2r]^ where 



N l (E , 2 ) = Re 



1 + e 2lVo -- 



1 - I «30,2 I 



1 - e 2 ^o, 2 j |2 _ e 2i(^o,2 | 2 ' 

are equal to the DOS, Eq. (|35| ) at energy So, 2- The resonance is slightly shifted from with 

5E lm) = J^L Im f 1 + ■ 1 + ^ 



(m) 



(B3) 



(B4) 



(B5) 



(B6) 



(B7) 



APPENDIX C: RESONANCE IN THREE-PARTICLE CURRENT 



In this appendix, the resonant form of the three-particle current, Eq. (pig), is derived. The M3o-matrix, which by 
definition is 



can be transformed using Eq. (Bl) to 

M 3Q = f-T- e i^f2ff-lf e ia, V tf-l 



D J WD 



(CI) 



(C2) 



21 



which can be written in the form 

M 30 = -4sin(^isin^ 2 (i/Vd+T^/D) + 2ia z sin (cpi + p 2 ) + e - ia * <p2 f- 1 e- i(T *' pl . 

\<f2 - mn\ <C 1 to 



It simplifies in the limit of D <C 1, \<f^\ = \fi — kn\ <C 1 and |^4"^l 



M 30 = 



{-l) k+r - 

L»3/2 



D - Ap[ k) p 2 m) ) (1 - a x ) + Dia z ( ^(1 - a x ) + ^(1 + a x ) 



,( fc ), 



(C3) 



(C4) 



Inserting this form of the M3o-matrix and the expansion (|33|) for r*o- and r 3 _|_ into Eq. ( j27| ) , as well as putting R = 1 , 
the probability current density for injection of a quasiparticle from the left takes the form 



[Mm 



(l-|a | 4 )(l-|a 3 | 4 )^ 3 

|l-e 2 ^o| 2 |l_ e 2iy 3 | 2 |Q| 2 



(C5) 



(m) l + e 2 ^ 



J _ g2i^ 



where Z? 1 is once again used. 

Since \<fi^\ <C 1 and |v4"^l ^ 1 anc ^ the DOS at energies -Ed, 3, Eq. (p5|), are equal to 



N (E) = Re 



1 + e 2itpo 



(C6) 



iV 1 ^) = Re 



1 + e 



2i¥>3 



J _ g2i(^ 3 



1-I«3| 4 

|1 - e 2 ^3| 5 



we arrive at the form 



[ME)] 1 = 



N l {E)N r (E 3 )D 3 



D. (fe) (m) 



(C7) 



(C8) 



APPENDIX D: RESONANCE IN FOUR-PARTICLE CURRENT 



In this appendix, we discuss the structure of the resonance in the four-particle current. The matrix 

M 40 = fe Ur '' P3 f- 1 e i "'' p *fe i ' 7 '' fil f- 1 , 

can be written as 



A/40 = 



D 2 



—8 sin tpi sin ip 2 sin p 3 yDT 1 + D sin (pi sin ip 2 sin 933 + D 2 sin (cpi + 933 — tp 2 ) 



(Dl) 
(D2) 



+2D sin c/?i cos ((^3 — p 2 ) VDT 1 + 2D sin</? 3 cos (<£>i — if 2) VDT 1 



1 

I) 2 



—AD sin t/?i sin 993 cos if 2 + 21? sin p 3 sin — (^2) V DT+ 



+2D sin <£>i sin (953 — ip 2 ) v DT +D cos ((pi + if 3 - p 2 ) 



From Eq. (D2) it is clear that in general M40 oc 1/D 2 . When both <pi and ip 3 are close to a multiple of n, M40 oc 1, 
while close to other double resonances, e.g. when ipi and ip 2 are close to a multiple of tt, M40 oc 1/D. 
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This results from the fact that the resonances are coupled with the MAR trajectory that crosses the barrier twice (see 
upper inset in Fig. |l^), instead of once as in the three-particle case, and therefore the resonance coupling is weaker. Another 
difference is that the width of the resonance and thus the height of the current peak is independent of the DOS at the 
entrance and exit energies, and only depends on the Andreev reflection probability at the exit and entrance energies. 
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